# This function computes the polynomial representing the syndrome from
# the received message bits
compSyndromePol := proc (rs, t, p, m, n, pol)
local S, pp, i, j, a;
    S := 0;
    a := RootOf(pol);
    for j to 2*t do pp := 0;
        for i from 0 to n-1 do pp := pp+rs[i+1]*a^(i*j) end do;
        S := `mod`(Expand(S+pp*x^(j-1)), p) end do;
    S
end proc;

p := compSyndromPol([1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1], 3, 2, 4, 15, alpha^4+alpha+1);

# The Euclid procedure returns the pair (v[i], r[i]) generated during
# the application of the standard Euclidean algorithm. The elements in
# the pair satisfy: degree(vi) <= mu, degree(ri) <= nu.
Euclid := proc (a, b, p, mu, nu)
local i, rs, qs, us, vs, ri, qi, ui, vi;
    rs := [a, b];
    vs := [0, 1];
    us := [1, 0];
    for i from 1 while not (rs[-1] = 0) do
        qi := Expand(quo(rs[i], rs[i+1], x)) mod p;
        ri := Expand(rem(rs[i], rs[i+1], x)) mod p;
        rs := [op(rs), ri];
        ui := Expand(us[i] - qi * us[i+1]) mod p;
        vi := Expand(vs[i] - qi * vs[i+1]) mod p;
        if (degree(vi) <= mu) and (degree(ri) <= nu) then
            break;
        end if;
        us := [op(us), ui];
        vs := [op(vs), vi];
    od;
    [vi, ri]
    # The following is the final result of the normal
    # Euclidean algorithm
    # [us[-2], vs[-2]]
end proc;

# Decoding function of the t-error-correcting BCH codes in GF(p^m)
# realized by polynomial pol.
# rs is the received message bits.
decode := proc (rs, t, p, m, pol)
local S, C, E, n, i, pr, v, d, a, ei;
    a := RootOf(pol);
    n := p^m-1;
    S := compSyndromePol(rs, t, p, m, n, pol);
    E := [];
    C := [];
    if S = 0 then
        print("No error occurs")
    else
        pr := Euclid(x^(2*t), S, p, t, t-1);
        v  := pr[1];
        d  := Expand(quo(v, eval(v, x=0), x)) mod p;
        for i from 0 to (n-1) do
            if ((Expand(eval(d, x = a^(-i))) mod p) = 0) then
                ei := 1;
            else
                ei := 0;
            end if;
            E := [op(E), ei];
            C := [op(C), (rs[i+1] + ei) mod p];
        od;
    end if;
    [C, E]
end proc;
